Langevin Equation for Particle in Thermal Photon Bath 
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The forward-backward path integral describing a charged particle moving in a thermal bath of 
photons is expressed in terms of the solution of a Langevin-type of equation. Approximate methods 
for solving this equation are discussed. 
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The interaction of charged nonrelativistic point particles with radiation is a fundamental problem of atomic physics. 
The electrons in light atoms move with nonrelativistic speed and may be described by a Schrodinger equation with 
the Hamilton operator 
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P- - A(x, t) 
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(1.1) 



where A(x, t) is the electromagnetic vector potential and V(x) = Za/r the Coulomb potential of the nucleus. The 
decay of an atom is governed by dipole radiation, and the matrix elements of the dipole operator between initial and 
final atomic states allow us to calculate directly the natural line width of a single atom in the vacuum. 

In actual physical systems, this simple situation becomes more involved. An atom may loose its energy in a variety 
of competing processes which lead to an additional broadening of its spectral lines. In this note we want to set up a 
theoretical framework for studying the broadening due to the interaction of the atom with a grand-canonical ensemble 
of photons in thermal equilibrium at high and moderately high temperatures. This will be referred to as a thermal 
bath of photons. Our work is a physically more realistic version of the well-known treatment of a particle in contact 
with a thermal bath of oscillators |l]-|| . 



II. DERIVATION OF THE LANGEVIN EQUATION 



The time evolution of the density matrix /9(x +a , x_ a ; t ) of the system described by (1.1) follows a quantum Liouville 
equation 

ihd tb p = H x+a p- H x _ a p. (2.1) 
This corresponds to a global time evolution equation 

p(x +a , x_ a ; t b ) = J dx +b dx- b U (x +b x_ b , < b |x +a x_ a , t a ) p (x +b , x_ b , t a ) . (2.2) 

The evolution kernel may be expressed as a forward-backward path integral due to Feynman and Vernon : 

U (x +b ,x_ b ,i b |x +a ,x_ a ,i Q ) = J Ox + Dx_i(x + (tt) - x +h ) 5(x-(t b ) - x_ b ) 

xexpj^' 6 ^(x^-x 2 _)-^(x + ) + nx-)-^x + A(x + ) + ^x_A(x_) J. (2.3) 

To describe a particle in a bath of photons, we perform a thermal average over the fluctuating vector potential A(x, t). 
This is done with the help of Wick's theorem, according to which the result can be expressed completely in terms of 
the time-ordered correlation functions at temperature T: 
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(Ti£ i (x,t)A J '(x',t / )) = G ij (x,x') = he 2 
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^ k r cosk(x — x') 



cos fik(i — t') coth — - — ^= — i sin flk(t — t') 



2k B T 



(2.4) 



where fik = c|k| are the frequencies of wave vector k, k B is the Boltzmann constant, and 5^ 



ij tr 



kikj/k 2 the 



transverse (5-function, We have used the four-vector notation x = (x, t) for space and time. The operator T is the 
time-ordering operator. 

At high temperature and for particle systems small in comparison to the mean wave length l/|k|, we may neglect 
the x-dependence in G{x,x'). Then G(x,x') can be approximated by 



G(x,x') 



2 h (k B T 
34^ 



5(t - 1') - 



\2{k B T) 



; 5(t - t') - 



i5(t - 1') 



(2.5) 



the brackets being the Fourier transform of the Taylor expansion of (hfl^ / 2k B T) coth hfl^/ 2k B T in (2.4). Moreover, 
the fluctuations around the average path x(t) = [x + (t) + x_(t)]/2 are small, so that we can approximate 



V(x + )-V(x_)~yW(x), 



(2.6) 



where y(t) = x+(t) — x_(t) is the difference between forward and backward paths. Then the path integral (2J3) takes 
the form (see j5| for more details) 



[/(x +fa ,x_ 6 ,i b |x +a ,x_ a ,i a ) = J VxDyd (x(t b ) - x b )5(y(t b ) - y b ) 



xex Pi - 



dt 
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y Mx - M7X + i— y - i 



wh 



2h J 2A{k B Tf 



yV^(x) 



(2.7) 



where we have abbreviated 



3nc 3 M' 



w = 2Mjk B T 



(2.8) 



At t a , the paths x(t) and y (t) s tart from x a = (x +a + x_ a )/2 and y a = x +a — x_ a , respectively. Representing the 
(5-function 8 (y(t b ) — y b ) in (2.7) as a Fourier integral, and inserting for y(t b ) the equation 



dt'y{t') + y a , 



(2.9) 



we rewrite the evolution kernel in the form 



U (x +b , x_ 6 , i & |x +a , x_ a , t a ) 
d3 P 1 n n 



(2tt) 3 

where we introduced the new variable 



dt 



iw 



y 6 Vl/(x) 



(yb -y a )t s ( x (^) _ x &) ■ 



T](t) = Afx(t) - Mjx(t) + J dt'VV (x(i')) - p. 



(2.10) 



(2.11) 



In Eq. ( 2.10| ), the variables y{t) at different points t a < t < t b are independent of each other, and we choose at the end 
points y{t a ) = y(t b ) = 0. This may be justified in a time sliced formulation, where the integrations over the variables 
next to the end points give only a trivial factor with respect to t he pro duct of integrals in which these variables are 
held fixed at the endpoint values. Furthermore, the path integral ( 2 . 1 d| ) does not depend on r)(t) outside the interval 
t € ( t a ,tb ), and is independent of rj(t b ) and rj(t a ). Hence we may choose r](t a ) — f](t b ) = 0, for convenience, and 
Eq. (2.11) can be solved as a differential equation on the interval t a < t < t b 



with the initial conditions 



Mx -M7X + W(x) = f](t), 
x(t a ) = Xa, Mx(t a ) - Mjx(t a ) = p. 



(2.12) 
(2.13) 
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We now perform the integral over y in Eq. ( fOOl ). We shall, from now on, neglect the expansion terms indicated by 
the dots, and obtain the path integral 



ta) 



d 3 p 



2?x exp 



1 f b 

* L dW{t) . 



ex Pi -- 



y b / dtVV(x) +p(y b - y a ) 



<5 (x(f&) - x 6 ) (2.14) 



where 77( f) de pends on x(t) via (|2 . 1 1| ) . There are some virtues of this representation in comparison with the path 
integral ( p. 10 ), in particular, if forward and backward paths start out and end at the same points, such that y a = 
y b = 0: the oscillatory integral in (2.10) is transformed into a Gaussian integral ( [2.14 ) which converges exponentially 
fast. Such a representation is obviously more suitable for numerical s imu lations. 

For convenience, let us express the evolution of the density matrix (2.2) in terms of the Wigner function defined by 
the Fourier transform 



W(x, p; f) = J d 3 y e^/ h p (x, y; t) 



(2.15) 



Here and in the sequel, we omit subscripts b from f for brevity. Then Eq. (2.14) yields the time evolution equation 
for the Wigner function as a functional integral 



W{x,p;t) = ^xexp ~J dt' V 2 (t') w(x(t),p- J dt' VF(x(t'));t 



(2.16) 



This equation has a simple physical interpretation. In the limit T — > 0, the functional integral ( 2.16] ) is concentrated 
arou nd r] (t') = 0, corresponding to a deterministic solution of Eq. ( [2.11 ) with T](t') = 0. In this limit, we obtain from 
Eq. (2.16) the Wigner function 



W{x, p; t) = W (x(f), p - J dt 1 VV{x{t'));t a 



(2.17) 



If the particle is decoupled from the bath, 7 = 0, we have p = Mx(t a ) and Mx(t) = p — f. dt VV r (x(t)), and we see 
that the time evolution of the Wigner function is given by the Liouville equation 



d t W + — p ■ V X W - V X V ■ V P W = 0. 



(2.18) 



The time evolution kernel can also b e expressed as a path integral over the noise variable rj(t). For this s imply 
change the integration variable in ( 2.16 ) from x(t) to 77(f). To find the functional determinant, we integrate Eq. ( 2.11 ) 
once more and write the result as 



dr 2 T7(f) = M[x(f) - x a ] - M 7 [x(f) - x(f a )] + J dsj ds' W(x(s')) - p(t - t a ). 



Differentiating this we obtain 

Sfii(t) 



Sxj(t r 



= Mdf 



M 



5ijS(t - t') - 7%<5(f - f') + — / ds / ds'\7 t \7 J V(x(s'))S(s' ~t') 



(2.19) 



(2.20) 



Applying the product formula Bet AB = Det ABetB, the identity Det (1 + K ) = exp (Tr K - ±Tr K 2 . . .) , and the 
prop erty Tr K n = 0, we see that the Jacobian for the transformation x(f) — > 77(f) is a constant, so that we can rewrite 
Eq. (|2T6[) as 



W (x, p;t) = (w (x(t), p - J dt 1 Vy(x(f')), fa 



(2.21) 



where the average with respect to 77(f) fluctuations is performed with the functional integral 



EE / DX 



exp 



~ I dtr, 2 (i) 
2w ' 



(2.22) 
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The calculation of (2.21) proceeds by solving first the Langevin equation (2.12) with the boundary conditions (2.13) 
to obtain the solution x(£), and subsequently take the expectation value with respect to the white noise with the 
correlation function 



{rf{t)ii>{t')) =wS ij S(t-t') 



(2.23) 



The quantum corrections to the Langevin equation are taken into account by replacing ( 2.23 ) by the colored-noise 
correlation function 



(?7 4 (£V(£')) = wF 3 coth 
= w5 lj 1 



ih d 



ih d 



2k B Tdt J 2k B Tdt 



\2{k B T) 2 dt 2 



S(t - t') 



S(t - t') 



(2.24) 



III. APPROXIMATIONS TO LIOUVILLE EQUATION 



From Eq. (2.23) we see that the average size of rj is \J2M^k B T and thus proportional to yT. Let us therefore set 



V — ■wT/ThV w hh the characteristic Bohr temperatu re Tj j = a 2 Mc 2 jk B (a = e 2 /Tic). Restricting ourselves now to 
one-dimensional systems, we search a solution of Eq. ( 2.12Q in the form 



(3.1) 



x (t) = [x cl (t) + v / T/T^Q(t)\ , 

where the first term on the right-hand side describes small dissipative correction to the classical equation 

Mx cl + V'(x cl ) = 0. (3.2) 



The second term on the right-hand side depends on the noise fj. Then, comparing terms of first order in 7 or ^T/Tji, 
and considering terms j^/T/Th as being of higher order, we obtain the following equations 



x cX B + 2Bx cX + B 



X c \ + — V" (X c i)x c i 



- jjV"(x cl )x cl = 0, 



(3.3) 



MQ + V"(x c i)Q = fj. (3.4) 

The solution x c \(t) depends on the boundary condition ( [?.13| ), in which we set 7 = in the present lowest-order 
approximation. 

The first equation is closely related to the second as can be seen by introducing A(t) = x c \(t)B(t) which satisfies 



MA + V"(x cl )A = V"{x cl )x cl . 



(3.5) 



Let G(t,t') be the Green function of the operator on the left-hand sides of (3.4) and (3.5), satisfying the harmonic 
differential equation 



dt 2 



n 2 (t) 



G(t,t') = 5(t-t r ), 



with the time-dependent frequency 



n 2 (t) = -v"( Xcl (t)). 



Then, the solutions of Eqs. (|3.3[) and (3.4) are given by 



Q{t)= / dt'G{t,t')fj{t') 
Jo 



(3.6) 



(3.7) 



(3.8) 



and 
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B V) = T^T* dfG(t,f)V"(x cl (f))x cl (t'), (3.9) 
Mx c i(t) Jo 

respectively. The Green function can be expressed in terms of two independent solutions £i(t) and ^(t) of the 
homogeneous version of the differential equation ( |Oi| ) dSIl] w ^ tn tne initial conditions £i(0) = and £i(0) = 1 and 
£2(0) = 1 and £i(0) = 0. These, in turn, can be obtained directly from the classical solution x c \(t) ||. 

If the classical solutions are bounded in all directions of phase space, then the Green function and its derivative 
are bounded. In such a case it follows from Eq. (3) in Ref. M that the random perturbation grows at most linearly 
in tT/Tu- A deviation of magnitude D from the classical solution appears at the time which is proportional to 
D/T. However, if the classical system is chaotic, which can happen either in multidimensional systems or in a one- 
dimensional system with a time-dependent potential, then the classical solutions as well as the Green function can 
grow exponentially fast in some directions of the phase space. In these directions, the classical solution expressed in 
terms of \JT /ThQ can grow to a size of the order D in a much shorter time ln(D/T). Such an increase of quantum 
corrections in chaotic systems has recently been discussed by H.Zurek p0| . For the quantum systems in a photon 
bath under discussion, this result is modified by the damping factor e~™\ Now, B(t) can be bounded, or it can 
increase linearly in t, as it happens for a linear oscillator. In a chaotic system, B{t) can also grow exponentially fast . 
Such a strong friction dampens completely the chaotic growth of quantum corrections, such that there is no need to 
apply quantum mec hanics to large macroscopic sys tems . It behaves like a strongly damped classical system. 

The solution (3.1) may be inserted into the Eq. ( 2.1 6| ) to find the time evolution of the Wigner function. 



Our semiclassical expansion can be extended systematically to any order, albeit with increasing complexity. Per- 
turbative as well as full nonperturbative solutions can be found most efficiently on a computer, yielding the full time 
evolution of the Wigner function and thus of the density matrix. 



IV. HARMONIC POTENTIAL 



The solution of our euqtions is simple for a harmonic potential 

Muj 2 



V(x) 



-x 2 . 



Then Eqs. (3.3) and Eq. (3.4) have the solutions 

B(t) = y (*-*«), <?(*) = ^ 



d£ smu){t-t')fi(t'), 



and (3.1) yields to lowest order in 7 the orbit 



1 • , IS 1 

— sin wis — t H 

1+jujMoj Mlo 



which determines the Wigner function via Eq. ( 2.17 ). 
The fluctuation width is 



(x(t)~(x(t))Y 



W f I 4-\- W 

—AG"*) = jm~2 



(t - i a )e^ 2 (*-*») 



dt' sinw(s - t')f)(t') 



sin 2uj(t — t a ) 



2w(i - t a ) 



(4.1) 



(4.2) 



(4.3) 



(4.4) 



For small times, this shows the same linear growth in time as for the Brownian motion of a free particle with u> = 0: 

((x(t)-(x(t))) 2 ) = -^(t-t a ), (4.5) 

As the time grows, the width oscillates around this behavior with frequency 2ui. For large times of the order 1/joj 2 , 
on the other hand, the width goes exponentially fast to zero (see Fig. |l|). 
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FIG. 1. Time dependence of fluctuation width in Eq. (4.4) for different values of friction constant 7. Time is measured in 
units of 1/uj. 



In the free-particle limit, the solution (3.1) reduces to 



(4.6) 



Inserting this into formula (2.16), we obtain the evolution of the Wigner function. For a wave packet with momentum 

(p—k) 2 a (x — x) 2 ' 



k and position x, the result is 



W(x,p;t a ) = -^exp 



Inserting (4.6) into ( |2. 16 ), we obtain for W the formula (4.7) with the replacements 



x — > x(t) = x — ~ to) an d a — > a(t) = a 



2w , 



(4.7) 



(4.8) 



i.e., a free evolution with the well-known spreading of the wave packet. 

For the oscillator potential, the mean position and the mean momentum of the particle run along the corresponding 
classical trajectories, and the thermal spread behaves as in (4.4). 

The spread in momentum space can be calculated from this using the relation 



((p(t)-«t))) 2 =MV 



(x(s) - (x(s))) 2 )ds. 



(4.9) 
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